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Abstract 

We apply the Lefschetz thimble formulation of field theories to a couple of different 
problems. We first address the solution of a complex 0-dimensional 0 4 theory. Although 
very simple, this toy-model makes us appreciate a few key issues of the method. In 
particular, we will solve the model by a correct accounting of all the thimbles giving 
a contribution to the partition function and we will discuss a number of algorithmic 
solutions to simulate this (simple) model. We will then move to a chiral random matrix 
(CRM) theory. This is a somehow more realistic setting, giving us once again the chance 
to tackle the same couple of fundamental questions: how many thimbles contribute to 
the solution? how can we make sure that we correctly sample configurations on the 
thimble? Since the exact result is known for the observable we study (a condensate), we 
can verify that, in the region of parameters we studied, only one thimble contributes 
and that the algorithmic solution that we set up works well, despite its very crude 
nature. The deviation of results from phase quenched ones highlights that in a certain 
region of parameter space there is a quite important sign problem. In view of this, the 
success of our thimble approach is quite a significant one. 


1 Introduction 


The so-called sign problem is one of the current big challenges for lattice field theories. It is 
in fact the major obstacle to tackling a non-perturbative study of the QCD phase diagram. 
Following pioneering work by Witten (1], Lefschetz thimble regularization has been proposed 
as a possible solution (2J [3] (for more recent contributions see also |U El El E|): the func¬ 
tional integral is defined in terms of fields taking values on non-trivial manifolds on which 
the imaginary part of the action stays piecewise (he. on the distinct thimbles attached to 
different critical points) constant. It is an elegant, although with many respects non-trivial 
alternative to the standard formulation of field theories. It has intriguing connections with 
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resurgence theory, a few results of which motivate the conjecture that the semi-classical ex¬ 
pansion of the path integral can be geometrized, as a sum over Lefschetz thimble 0 • Morse 
theory [9] is the natural framework for discussing the thimble regularization, even though it 
could be that it does not necessarily have the last word on the subject. In this work we will 
discuss the thimble solution of two different models, having in mind two big issues. First of 
all, since there is a thimble attached to every critical point of the (complexified) theory one is 
considering, we need to understand how many thimbles do give contribution to the solution 
of a given theory. A second relevant problem is that of devising Monte Carlo algorithms to 
correctly sample thimbles, which are manifolds for which we lack a local description. 


Toy models can be a precious tool to approach hard problems in theoretical physics, with 
the hope that a simplified model can nevertheless capture the relevant issues. Being the sign 
problem both relevant and hard, it does not come as a surprise that toy models have been 
around for a while. Quite interestingly, some of these have been resisting the efforts to solve 
them in much the same way the real problem is still far for being fully solved. The first ap¬ 
plication we discuss is the solution of a toy model that dates back to almost thirty years ago 
[TO] . It can be regarded as a 0-dim version of a <f A held theory. This simple model became a 
benchmark for the complex Langevin treatment of the sign problem, and quite interestingly 
only partial success has been claimed over the years. We will show that the thimble regu¬ 
larization completely solves the model. In this case we will show how in different regimes 
one or more thimbles give contribution to the solution. Moreover, it is possible to perform 
numerical simulations on the thimble(s): in this simple case we will have a number of al¬ 
gorithmic solutions at work. A partial account of these results has already been given in HU- 


We will then address the solution of a chiral random matrix theory. This is a somehow 
more realistic problem, for which once again the application of the complex Langevin method 
has been shown to be non-trivial: in a given parametrization it fails [12], while in a different 
one it works Actually one can show that the sign problem can be quite severe for this 
model. The theory has an adjustable parameter (the dimension N of the matrices) which 
controls the dimensionality of the problem one has to solve. Since the analytical solution 
for the observable we will study (a mass condensate) is known, it appeared to us a perfect 
setting for testing a conjecture: it could well be that more than one thimble contribute (just 
like in the 0-dim 0 4 theory) in low dimensions, but as N grows it could be that a single 
thimble dominates in the thermodynamic limit. While we were ready for a richer scenario, 
in the region of parameters we studied we actually did not find any other thimble but the one 
attached to the global minimum. We did not find any problem related to the parametriza¬ 
tion of the theory; in particular, the parametrization that was failing in the case of complex 
Langevin was absolutely fine for the thimble treatment of the theory. In this case algorithmic 
problems are non-trivial. We will show how a natural parametrization of the thimble can be 
the starting point for an algorithmic solution that in this case works in its simplest version 
(admittedly a very crude one). 


The paper is organized as follows. In section 2 we give a brief account of the Lefschetz 
lr This is a literal citation from [5]; we will have more to quote later on the subject of resurgence. 
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thimble approach to field theories: this is a short review of results that are collected to fa¬ 
cilitate the reader. Section 3 is dedicated to the 0-dim </> 4 toy model, showing that thimbles 
provide a complete solution; in particular we show that we can effectively numerically simu¬ 
late the model, making use of different algorithms. In section 4 we address the chiral random 
matrix theory, showing that we can numerically solve it by thimble regularization: all the 
analytical results are correctly reconstructed. In the final section we draw a few conclusions 
and mention natural steps forward. 


2 Thimble regularization in a nutshell 


In the following we collect the basics of thimble regularization: the interested reader is re¬ 
ferred to wot for further details and references. 


A conceptual starting point to approach the thimble regularization is that of generalizing 
saddle point integration. The latter displays a couple of features which appear as good can¬ 
didates to tackle the sign problem: stationary phase and localization. The full generalization 
of saddle point techniques is formulated in the framework of Morse theory |3j. 

• One starts with an integral on a real domain of the form 

d n x g(xi ,..., x n ) e _5 0ir--^n) (1) 

in which C is a real domain of real dimension /^] and both S(x i,... , x n ) = S R {x i, ..., x n )+ 
iSiix i,... ,x n ) and g(x i,... ,x n ) are holomorphic functions. The notation for the ex¬ 
ponential makes it clear that we have already in mind a functional integral (even if the 
normalizing factor Z~ L is missing). For such an integral the following decomposition 
holds 

[ d n xg(x 1 ,...,x n )e- s ^’-’ x ^ = Vrv [ d n zg(z 1 ,...,z n )e~ s ^’-^ (2) 

JC a Jj a 

in which an extension from a real domain to a complex one has been performed (see 
the complex variables z t — x { + ii/i as opposed to the real ones Xi). (j2]) holds in the 
homological sense, i.e. C = '^ Ja n a J a . 

• The index a labels the critical points of the complex(ified) S(z\,... ,z n ) and to each 
critical point p a a stable thimble J a is attached. Each J a is defined as the union of all 
the Steepest Ascent (in the following SA; we will also write SD for Steepest Descent) 
paths falling into p a at (minus) infinite time, i.e. the union of the solutions of 

dxi 
dr 
dyj 
dr 

2 In what follows we will be a little bit sloppy in our notation: whenever there is not subscript attached 
to a symbol (e.g. x), that will denote a n-dimensional coordinate. 


dS R (x,y) 

dxi 

dS R (x,y) 

dyi 


(3) 
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satisfying z(r = — oo) = x{r = —oo) + iy(r = —oo) = p a . The real dimension of each 
thimble J a is n. It is quite natural to regard thimbles as manifolds embedded in C n 
(which is instead of real dimension 2 n). 

• For each critical point one also defines unstable thimbles KL a as the union of all flows sat¬ 
isfying equation ([3]) and going to the critical point in the opposite time limit, i.e. such 
that z(r = oo) = p a . The coefficients n a count the intersections of the /C CT with the 
original domain of integration n a = (C,/C CT ). 

• The imaginary part Sj stays constant on thimbles, i.e. there is a phase associated to 
each thimble. 


Note that in the framework of field theories a natural picture of universality emerges. A 
single thimble can give us a formulation of a field theory with the same degrees of freedom, 
the same symmetric^] and symmetry representations, the same Perturbation Theory and 
naive continuum limit of the original formulation (see [2] for details). In force of universality 
we expect that these properties essentially determine the behavior of physical quantities in 
the continuum limit. Moreover a simple argument suggests that in the thermodynamic limit 
only thimbles attached to global minima can survive, as it is easily seen (we now call (f a the 
critical points, having in mind field configurations, and consider the partition function of the 
field theory) 



In the end, it could well be that a full resolution in terms of all the thimbles could turn out 
to be with many respect overkilling the original problem. Having said this, we nevertheless 
stress that both the universality argument and the thermodynamic argument can not be re¬ 
garded as conclusive. It is worth noting that resurgence theory m with many respect even 
asks for more than one thimble in view of the interpretation of the semi-classical approxi¬ 
mations as trans-serie^} All in all, the fact that in certain cases a single thimble dominance 
can take place has to be regarded as a conjecture: as we will see, this was in a sense one 
of the motivation of this work on a Random Matrix Model. The subject deserves deeper 
investigation and we think it will certainly receive it. 


It is good to have a somehow more constructive approach to the thimble formulation. 
We therefore now sketch a few more technical details that the reader will see at work in 
the following sections. The integral we have in mind will be the functional integral of a 
field theory. Let us first of all parametrize the field in the vicinity of a critical point as 
= 4>i — Here and in the following i is a multi-index; in particular it can refer to a 

3 Since we have discussed the case of a non-degenerate hessian, one could wonder how the method can 
be applied in case where Spontaneous Symmetry Breaking (SSB) is in place. In |TJ[ a solution has been 
described and shown to be effective: one introduces an explicit Symmetry Breaking term h and studies the 
limit h — » 0. This is not the only way to proceed: for a different thimble approach to SSB the reader is 
referred to [3j. Symmetries are dealt in yet another way in the case of Gauge Theories: the construction of 
thimbles was discussed in [15j and reviewed in [T] ; see also the discussion in mm- 

4 A nice account of many issues connected to resurgence has been recently provided in [18] , 
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real or imaginary part. The real part of the action can be expressed as 


Sr (0) = S R (0 ff ) + -<L r iK> + O (0 3 ) 


where the 2 n x 2 n matrix H is the hessian evaluated at the critical point 

d 2 S R 




4* —^(7 


( 4 ) 


H can be put in diagonal form 


A = diag (Ai, • • • , A n , -Ai, • • • , -A„) 

by a transformation H = WAW T defined by the orthogonal matrix W whose columns are 
given by the normalized eigenvectors of H , that is {r/4} i=1 ... 2n . Half of the eigenvalues of 
H are positive; the corresponding eigenvectors span the tangent space to the thimble at the 
critical point. Any combination of these vectors is a direction along which the real part of 
the action grows. If we leave the critical point along these directions integrating the SA 
equations we span the thimble. On the other side, the other directions (which are attached 
to the negative eigenvalues) would take us along the unstable thimble. 

At a generic point Z £ J a we miss a priori the knowledge of the tangent space TzJa ; hr 
general we expect that the latter is not parallel to the canonical basis of C n whose duals 
appear in d n z = dzi A • • • A d z n . We thus want to perform the relevant change of coordi¬ 
nates from the canonical ones (of C n ) to the basis of TzJa , given by the (complex) vectors 
{U^}i = i„. n (these are orthonormal with respect to the standard hermitian metric of C n ). 
Let ip : N C J a — > [R n be a local chart in a neighbourhood N (Z J a oi Z 




z+j2 u(l 


i =1 


Y + O ( y 2 ) e 


If we denote U the n x n complex unitary matrix whose columns are the vectors {U^}, we 
can express the integral of a generic function / (Z) on the thimble as 


J d"zf(Z) = I deter p.-pr)) 

N ‘f(N) 1=1 


(5) 


In this expression the quantity det U = e tuJ (U is unitary) has appeared; this is what has 
been termed the residual phase [2J (see ra for further details). This could in principle re¬ 
introduce a sign problem in the thimble formulation, but it is expected that this is not the 
case. Not any phase gives raise to a serious sign problem, and in particular one expects that 
a phase changing rather smoothly can be safely taken into account by reweighting. This 
expectation could appear optimistic, but has been till now confirmed (see 0 ) and will be 
confirmed also in this work. 


We end this brief introduction to the thimble formulation by going back to the construc¬ 
tive point of view: we can span the thimble by integrating the SA equations for the field 


5 




^0 


i — 1 ■ ■■2n 


d0, _ dS R 
d t d(j)i 

This has a counterpart in parallel-transport equations for the n basis vectors which defines 
the tangent space to the thimble (see BED 


dVf = - w 

dt “ dcpkdcpj k 


i = 1 • • • n j = 1 • • • 2 n. 


( 6 ) 


We can set up a similar equation for any other vector with an initial condition on the tangent 
space at the critical point; © expresses the parallel transport of a vector along the gradient 
flow. In the vicinity of the critical point one knows the asymptotic (t —> —oo) solutions 


t < 1 


n 

(f)j ( t ) « (f> a j + Y Vj e Xit rii j — 1 • • • 2n \n\ 2 = 1 

i= 1 

f/dd (t) ~ v^e Xit j = 1 ■■ - 2n i — 1 • ■ ■ n 


(7) 


Note that from a practical point of view the former parametrization is viable only provided 
one introduces a reference time to 1 at which the former asymptotic solution holds. We 
will make estensive use of the former equations, which in particular can be regarded as initial 
conditions for a given flow on the thimble, e.g. for eq. 0- It thus emerges a natural picture 
in which a generic point T £ 3<j is unambiguously defined by a choice of €I and the time t : 
$ = <3> (h,t) (this has been very effectively discussed in [3j). Note also that one could insist 
on regarding the ([7]) as valid all over; this would in turn mean one is considering a purely 
quadratic action (he. the free field approximation). 


3 The 0-dim 0 4 toy model 


We will now put at work what we have just seen, applying the thimble regularization to the 
study of the action 

S (0) = + ^A0 4 

with 0 G [R, A G [R + and a = cr R + z ay G C. This is obviously a toy model, and one regards 
as correlators plain one-dimensional integrals such as 

(0”) = t/#^"e- s W (8) 

R 

with the partition function given by: 

Z = j d0e- swo 

R 

5 It is worth to recall here that the subscript i is a multi-index, in which real and imaginary parts are on 
the same footing. 

6 The hat notation remembers us of the normalization condition, i.e. fi singles out a direction in the 
tangent space at the critical point. 
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The solution is given in terms of a modified Bessel function, i.e. Z = yA^esx ^(^’dif¬ 
ferentiating appropriately which one can get any of the (|8|. The choice of a complex a is a 
prototypal case of the sign problem: with a complex action, we miss a positive semi-definite 
measure and hence a probability distribution to start with; in particular, a direct access to 
Monte Carlo methods is ruled out. 


It was realized long time ago that a solution to the sign problem could be searched in 
the context of Stochastic Quantization: the Langevin equation admits a formal solution 
also for complex actions, in particular via the Fokker-Planck formulation |2U| EL] • Turning 
the formal arguments into a rigorous proof eventually turned out to be hard and numerical 
instabilities (suggesting problems) were in particular discussed in the context of the theory 
at hand ra. Much experience has been gained over the years and much progress has been 
done [223. The question of convergence of complex Langevin equation remains a subtle one, 
and quite interestingly even the simple model at hand displays delicate issues. For a recent 
and thorough study of complex Langevin dynamics of this model, the reader can refer to 
[23]. One peculiar feature of this model is that complex Langevin simulations display diver¬ 
gences for (c i > n ) with n > 4 in a certain region of parameters. The relation between complex 
Langevin and thimbles has been investigated in [231125]. 


We can complexify the field by setting 0 — x + iy. As a result, real and imaginary part 
of the action read 

S R = ^ [a R (x 2 -y 2 ) -2(nxy] + ^X (x 4 + y 4 - 6x 2 y 2 ) 

s 1 = \ [ a i ( x 2 - y 2 ) + 2 a R x v\ + A ( x3 y - x y 3 ) 

The hessian is built from the second derivatives of Sr and takes the form: 


H(x,y) 


f&R + 3Xx 2 — 3A y 2 ■—(T x — 6A xy A 

y — a 1 — 6A xy —(Jr — SXx 2 + 3A y 2 J 


(9) 


There are 3 critical points: 0o = 0 and (j)± = (which are the two, complex valued, 

“Higgs vacua”). The question is now which thimbles do give a contribution to the integrals 
we want to compute, and the answer is quite different in the 3 cases a R > 0, or < 0 and 
ctr = 0: in each case we computed the stable and unstable thimbles associated to each 
critical point. This can be done putting at work the constructive definition of thimbles we 
discussed in the previous section. 

I 11 practice, we want to integrate the equations of SA starting in the vicinity of the critical 
point for an arbitrarily long flow time t. We can do this provided that the initial condition 
is chosen correctly: for the stable thimble this means we leave the critical point along the 
direction (in the xy plane) which is given by the eigenvector of positive eigenvalue of the 
hessian ([ 9 ]) computed at the critical point. Once we have singled out the relevant direction, 
we can ascend in two ways (namely, increasing or decreasing x), both of which we have to 
take to cover all the thimble. By holomorphicity the hessian has two eigenvalues opposite 
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in sign. Since Sr always increases along the flow, exp (—Sr) goes to 0 as t —» +oo, thus 
ensuring convergence of the integrals along the thimble. To obtain the unstable thimble /C^, 
we can repeat the same procedure described above, but picking up the eigenvector of the 
hessian of Sr with negative eigenvalue. Note that the unstable thimble is needed because 
the coefficient n a in our master equation (J2| counts the intersection of such thimbles with 
the original domain of integration, which in our case is the real axis (the sign ambiguity is 
not resolved just by this definition, but it can be deduced by means of other considerations). 
Figures [I] (left panel) and [2] show the results for the three cases <j r > 0, ctr < 0 and or = 0 
(see also [2B]). 




Figure 1: Thimbles structure for a = 0.5 + zO.75, A = 2 (left panel). In this case only 
the unstable thimble attached to z = 0 intersects the real axis and thus only one critical 
point contributes. On the right we can see how the Langevin simulation correctly covers the 
relevant thimble. 

From figure [lj we see that when or > 0 the unstable thimbles related to the Higgs vacua 
do not intersect the real axis. Therefore these points do not contribute to the integrals, that 
is n± = 0 and Uq = 1. By integrating along the stable thimble attached to (ft o, we recover the 
correct results for, say, Z = f e~ s (the integration can be easily carried on along the real axis 
both analytically and numerically). The case (Jr < 0 depicted in the left panel of figure [2] is 
a totally different matter, as we cross the Stokes ray ctr = 0 while changing sign to ctr. Now 
we see that the unstable thimbles connected to the Higgs vacua do intersect the real axis and 
therefore n± ^ 0, as well as no ^ 0. The correct combination which recovers the expected 
results for the integrals turns out to be n 0 = —1 and n± — +1. What is the origin of this 
discontinuity? and, above all, if we hadn’t known the correct result from the beginning, 
how would have we calculated the n a l The answer lies in considering the case ctr = 0, 
showed in the right panel of figure [2} The stable thimble connected to 0 exhibits the Stokes 
phenomenon: in fact it “collapses” into the Higgs vacua, from which it does not “move” any 
more; the unstable thimble continues to say that n 0 ^ 0. The stable thimbles connected 
to the Higgs vacua display the same shape, but their unstable counterparts collapse into 0 
(by overlapping its stable thimble) and therefore there is intersection with the real axis; so, 
n± 7 ^ 0. However, there is no integer-valued combination of n a that recovers the correct 
results for ctr = 0. This is quite expected, as the Morse decomposition along thimbles is 




































































not legitimate when we are on a Stokes ray, on which we clearly are (the imaginary axis in 
the complex a plane is a “Stokes ray” Now, the original integral is continuous (in fact, 
it is holomorphic) in a and therefore there cannot be any discontinuity in the computation 
of the partition function Z in a R = 0. Thus, we must have Z [ctr —» 0 + ] = Z [a R —> 0 _ ] = 
Z [ctr = 0]. By examining the integration along the thimble connected to 0, we find that it 
is discontinuous in a R = 0, and again, this is not surprising as the thimble shape undergoes 
a radical change between the two cases. The change in sign of the n a is precisely the only 
one which keeps the original integral continuous while crossing the Stokes ray. 




Figure 2: Thimbles structure for a = —0.5 + i0.75, (left panel) and a = 70.75 (right); in 
both cases A = 2. For a R < 0 (left) all the three critical points contribute. a R = 0 (right) is 
an example of a Stokes phenomenon. 


3.1 A variety of algorithmic solutions 

Within the thimble regularization we were able to perform numerical simulations of the quar- 
tic toy model, making use of different algorithms. In particular, we were able to numerical 
compute all the possible moments ([8]). 


It was observed in (2j that the Langevin algorithm is the obvious candidate for sampling 
configurations on the thimble. In 


d fc _ dS R 

~dF~~ay + ,,i 


i — 1 • •• 2n 


( 10 ) 


the drift term constrains the held on the thimble by definition, so that the problem boils 
down to extracting a convenient noise, i.e. a noise tangent to the thimble. We do not discuss 
here the original solution which was put forward in [2] (the Aurora algorithm); there will 
be a convenient time for such a discussion when we later approach the CRM model. Here 

7 see T] for a detailed explanation of the Stokes phenomenon with respect to the Airy integral 
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it suffices to say that, being the thimble 1-dimensional, at every point the tangent space 
reduces to the direction singled out by the drift term itself. As a matter of fact, Langevin 
works pretty well; in the right panel of figure [l] one can see how the simulation correctly 
samples configurations on the thimble. Here parameters are the same of the left panel, so 
one thimble is relevant, i.e. the one attached to the origin (which in the notation of (|2]) 
we denote po). The (Aurora) Langevin algorithm samples points according to the measure 
normalized byj^] 

Z (0) = [ dr e~ SR (11) 

■ho 


We now denote 


(O)o = 


f Ja dr O e Sr 
ZW 


( 12 ) 


and stress that this is not what we have to compute. Properly including the residual phase, 
the correct result was computed as 


(O) 


(e^O) o 

(e^)o 


(13) 


When < 0 the thimbles associated to all the three critical point^] contribute and we 
have to compute 


(O) 


Y^i=o n i e 1 S,( ' Pi ' > fj d re Sr O e li 
Y^=o n i e~’ lSl ( pi ') Jj dr e~ SR e luJ 


(14) 


which can be written 


_ (e iu O )o + ai(e iui O)i + a 2 {e iuj O ) 2 
' ~ (e iuj )o + a 1 (e i “) 1 + a 2 (e i “)2 


(15) 


with 


n . e -iS 7 ( Pi ) Z p) 

ot{ = - . g , ' i — 1,2 (16) 

n 0 e- iS hpo)2(0) ’ v ’ 

On each thimble J % (i = 0,1, 2) the quantities ( e lul O)i and ( e lu )i can be computed via 
(Aurora) Langevin simulations. The (complex) unknown coefficients cq can then be fixed by 
relations which can be regarded as renormalization conditions in a physical scheme, i.e. 


{e lu> Oj) o + ai(e itxJ Oi)i + a 2 {e luJ Oi) 2 
(e iu} )o + oti(e iu )i + a 2 (e iu} ) 2 

8 r is the real coordinate on the (one-dimensional) thimble. 
9 In the notation of |2| we now denote <f> + = p\ and (j>- = p 2 - 


(17) 
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where the X t are known values of given observables Oi ( e.g. , in the case of moments 
(|8j) , two of them). As always in such an approach, one gives up predicting everything, but 
after normalizing results to a (minimum) number of external inputs, one has full predictive 
power for (all the) other quantities. Of course computing the moments ([8]) for the toy model 
at hand is not such a big numerical success; nevertheless the outline of the method is quite 
general. In particular, we will refer to it in section 4.3. 


Another algorithmic solution for this simple setting is provided by the Metropolis algo¬ 
rithm which is described in [27] , The method relies on a correspondence between the full 
model one has to simulate and a gaussian approximation associated to it. The latter is 
obtained by diagonalizing the hessian at a critical point and truncating the expansion of the 
action around it, i.e. 


2 n 


SrW = Sn(M + R Y1 + s a(’l) 


(18) 


k=1 


where the rjk are the = 4>k — 4>a,k of equation Q expressed in the basis provided by the 
eigenvectors of the hessian, with a convenient ordering in which A& > 0 for k — 1... n and 
A*, < 0 for k = n + 1... 2 n. For the gaussian action (18) it is very simple to construct the 


associated stable thimble. It is a flat thimble in which the tangent space is known once and 
for all, i.e. the span of the eigenvectors associated to {A ^|k = 1... n}: we term it a gaussian 
thimble. For the gaussian thimble the solution in the right hand side of ([T]) is valid all over 
the manifold. 


The simulation is run as a quite standard Metropolis algorithm controlled by an ac¬ 
cept/reject test, with a mechanism for proposing configurations which is dictated by the 
correspondence between the thimble one has to sample and its gaussian approximation. We 
sketch the method in the case of more than one thimble contributing to the final result, to 
stress how also in this case we were able to run numerical simulations on thimbles, for both 
<7 r > 0 and < 0. 

The method always handles a couple of configurations, i.e. one (j) field on the thimble 
we have to sample and one auxiliary r] field on the associated gaussian thimble. In order to 
extract a new field one proceeds as follows: 

• One proposes a thimble a' (i.e. a critical point) with a probability 

\n a '\ 

Ea Kl 

• One extracts a configuration rf on the gaussian thimble associated to that critical point 
according to the weight e~ s °. This is trivial, given the gaussian form. 

• One starts a SD on the gaussian thimble with 7 / as initial condition. The integration 
is carried on over a time extent t such that one ends up close enough to the critical 
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point, namely at a point where the gaussian thimble and the thimble one has to sample 
effectively sit on top of each other (this means that the solution ([7]) holds for both 
thimbles). We call 77 the configuration that has been obtained in this way. 

• Taking 77 as the initial condition, one integrates the SA equations for the complete 
theory over the same time extent f. This generates the new configuration (j)'. 

• 0 ' is accepted with probability 

P acc = min jl, e -[ s « ( ^')- 5 i i W]+[ 5 G (r / ')-S G ( 7 ? )]| 


The result for a given observable O is obtained as 

I e‘“<« 


where the index t runs over all the configurations sampled by Metropolis. 


Notice that the previous accounting of the Metropolis algorithm is technically different 
from the proposal of [27]p^j The latter relies on an exponential mapping for the integration 
time (he. r = e - *) and an adjustable parameter is introduced to control convergence prop¬ 
erties (the interested reader can refer to figure 4 of [IT]). For a given, effective choice of this 
parameter figure [3] displays how the three thimbles giving contribution in the region (Jr < 0 
are sampled in a Metropolis simulation. 


We stress that also in this case one could think of situations in which the weights n a are 
unknown. However, here the situation is different from that of Langevin, since we only need 
to know a few integers values. I 11 other terms, given the knowledge of the set of relevant 
integers {n a }, the problem is solved on an entire region of parameter space: in the case at 
hand, for each ctr < 0 (technically, over the entire region which ends up in a point where a 
Stokes phenomenon shows up). Notice that in principle there can be different ways of finding 
the relevant set of integers ( e.g. known asymptotic solutions in a convenient region) ^ 


Both (Aurora) Langevin and Metropolis could correctly compute the moments (J 8 |, in 
both regions Or > 0 and <Jr < 0. For example, figure [3] (right panel) displays the computed 
values of ( 0 8 ) over a range of both ctr > 0 (one thimble being relevant) and ctr < 0 (three 
thimbles to be taken into account) values. 


Note that there is another natural way of computing on a thimble, and this takes ad¬ 
vantage of the fact that on a thimble there is a one-to-one correspondence in between con¬ 
figurations and values of Sr. To make things simple, let us consider the case in which only 
one thimble is relevant and let us write the partition function Z = Z up + Z down : these are 

10 We decided to enlighten the rationale of the algorithm, leaving out the technicalities. 

11 The cn introduced for the (Aurora) Langevin algorithm entail instead the values of partition functions 
and are given at a given point of parameter space; in the case at hand, for a given value of ur < 0. 
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Figure 3: Left panel: the three thimbles associated to a = —0.5 + 10.75 correctly sampled 
by a Metropolis simulation. Right panel: for the same choice of parameters, the computed 
values of (</> 8 ) over a range of both or > 0 and or < 0. 


the two contributions resulting from the two pieces of the thimble we have already referred 
to. Namely, they are associated to leaving the critical point along the direction dictated by 
the eigenvector of the hessian in one of the two possible ways (i.e. increasing or decreasing 
x values). Each Z up/down has the global phase e~ lSl(j> ^ as a factor and features an integrand 
which is the residual phase times a monotonic function of Sr. It thus can be written taking 
the action as the integration variable, e.g. 


Z = e-iSfr,) 


dS 


,-s lt 


VS 


R 


-1 i tan 1 (d y S R /d x S R ) 


(19) 


We could have written the integral by taking the flow time as the integration variable (also 
in this case there is a one-to-one correspondence with the configurations along the thimble, 
each reached at a given flow time). Note that in computing (19) one proceeds by integrating 
the SA. We illustrated the issue by taking into account the Z, but we showed that all the 
moments (J8]) can be successfully computed in this way. In a sense, (19) is the prototype of 
a parametrization we will see at work for the CRM model. 


All in all, we think that the simple toy model we discussed is a perfect playground to see 
thimble regularization at work: it is instructive both from the point of view of inspecting 
the structure of relevant thimbles and from the algorithmic point of view (we can compute 
on thimbles). 
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4 Chiral random matrix model 


We now address the chiral random matrix model defined by the partition function 



J did't det Nf (D{p) + m) exp (-N ■ Tr[4 , 4 + 4*4]), 


( 20 ) 


where 


.D(/i) + m 


m 

i cosh(/i)$' 1 ' + sinh(/i)'ld 


i cosh(/r)<f) + sinh(/i)'P 
m 


( 21 ) 


The degrees of freedom of the model are N x N general complex matrices T and <f>. Since its 
introduction it has attracted attention due to the many features which it shares with QCD 
[28 , 29j [30]: they both have in their functional integral the determinant of a Dirac operator 
and the flavor symmetries and explicit breaking hereof are identical. Chiral perturbation 
theory at leading order in the e-domain is the relevant low energy theory in the microscopic 
limit for both theories, which resulted in a lot of interesting insights into QCD coming from 
the (much simpler) random matrix theory. The microscopic limit in which contact is made 
with e-regime of chiral perturbation theory is that of IV —)■ oo with m = Nm and fi = y/N /i 
kept constant. 

A sign problem is there for this theory as it is for QCD. This sign problem can be a severe 
one, as it is made manifest by considering the observable we will be concerned with, i.e. the 
mass dependent chiral condensate 




i°g (Z) 


( 22 ) 


Figure [4] displays both the exact (solid red line) and the phase quenched (dashed blue line) 
results for A ( 7777 ) as a function of fh for N = 1,2,3,4 and fixed Nf — 2, fi — 2. As it can 
be seen, the sign problem is indeed severe in certain regimes of small (rescaled) masses 1 ^ 


Our interest in the model was triggered by [ 12 , IT3] : the nature of the sign problem 
(which is due to the determinant) has a counterpart in a non-trivial success of the complex 
Langevin method (which needs to take the logarithm of the determinant to define the stan¬ 
dard effective action dictating the drift term of Langevin equation). While the application 
of complex Langevin in the most direct parametrization of the theory fails ra. a different 
parametrization (resulting in a different complexihcation) reproduces the right results [ 13] . 

4.1 How many thimbles should we take into account? 

We take the most direct path to complexification, i.e. for each field (we directly deal with the 
matrix elements) = a %] +ibij and = a ij +if3ij, each real component gets complexified 

( e.g. fiij = + ifijp)- We adhere to the notation of [T2] and denote the action as 

S (a, b,a,P) = Njp «• + bl + a% + /%) - N f Tr log (m 2 l NxN - XY) 

12 The value of the condensate is real. This has to be understood later when we will compare to our results; 
we will always plot only the real part, the imaginary one having been correctly verified to be zero within 
errors. 
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Figure 4: Exact (solid red line) and phase quenched (dashed blue line) results for the con¬ 
densate, at fixed Nf — 2, jl — 2. 


with 


Xij = i cosh fi ( a,ij + ibij ) + sinli fi (a^- + ify) 

Yij = i cosh fi ( a,ji — ibji ) + sinli fi (<x,j — i/3ji) 

Once we have complexified the degrees of freedom, the Erst step for the thimble approach 
is the identification of critical points of the resulting action. First candidate is the absolute 
minimum which is already there for the real formulation, i.e. T = <3> = 0. All the relevant 

formulae for the spectral analysis of the hessian of Sr are collected in appendix [Aj Here we 

simply state that the hessian in 0 has the expected number of positive eigenvalues, i.e. the 
real dimension of the thimble attached to 0 is 47V 2 . Note that there is a huge degeneracy: 
we have only two different eigenvalues, with the two eigenspaces having the same dimension. 
As the (rescaled) mass fh gets smaller, the gap between the two eigenvalues gets larger. 
Some insight can now be gained from equation ([7]): in first approximation, the closer the 
eigenvalues, the more isotropic we expect the thimble to be. This expectation turned out to 
be correct in view of the results of our simulations. 


We tried to identify other critical points. In our study we explored different values of fh 
(at different values of N) while keeping fixed Nf — 2 and jl — 2. One approach was solving 
VS = 0 via Newton-Raphson method. We cross-checked results by applying the Nelder- 
Mead simplex method to minimize || VS || 2 . We found two classes of extrema, both outside 
the original domain and featuring an action smaller than Sr(^ = $ = 0), which turns 
out to be the absolute minimum in the original domain. Under such conditions, since the 
unstable thimbles attached to the extrema we found can not intersect the original domain of 
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integration, we expect no contribution from their stable thimbles (n a = 0; see section II.B.3 
of pj for a more extensive discussion). 


4.2 Algorithmic issues for the CRM model 


While for the 0-dim toy model the original algorithmic solution proposed in [2] is trivial, this 
is not the case for the CRM model. However, previous experience with the Bose gas [Til 
taught us that there can be lucky cases. Let us remind the reader of the Aurora algorithm 
and of its gaussian approximation (which successfully deals with the lucky cases we were 
referring to). 

We want to extract a proper noise vector for the Langevin dynamics 


d (j)i _ dS R 

n ‘ 


i — 1 • •• 2n 


We can proceed as follows [2]: 


• We extract a gaussian noise (0), where the superscript qualifies this quantity as an 
initial proposal and the argument has to be thought as a flow time in a sense that will 
be clear soon. 


• We evolve it following the flow (|6]) downwards (i.e. with a change of sign with respect 
to (JbJ), aiming at getting close enough to the critical point in order to make contact 
with the regime of ([7]). This will hold at a given descent time r*. 

• We then project with 


Vi 





and normalize the result 



r being extracted according to the n-dimensional x distribution. 


(23) 


(24) 


• We then ascend along the flow, covering again a time interval of length r*. The result 
is the noise rj, we will put in our Langevin equation. 


Extracting the noise vector is not yet the end of the story, since any finite order approximation 
to Langevin equation, e.g. the Euler scheme 

! _ 3S}i rz~ 

4>i = 4>i~ St — + Vtit rji 

will introduce systematic effects; since the manifold is not flat, the final point 0' will be 
moved away from the thimble. The obvious remedy for this effect is to repeat just the same 
we did for the noise vector (move the configuration along the flow downward, close to the 
critical point, project it onto the tangent space, move it upward along the flow for the same 
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time length). Note that it is expected that, in all the descent/ascent mechanisms we have 
just described, the downward flow, i.e. the SD will be numerically delicate. It is thus much 
better to formulate the descent as a boundary value problem (BVP) rather than as an initial 
value problem, as it was observed in [IT] , 


A much more appealing observation was also made in m there are lucky cases in which 
a quite rough approximation holds; with a slight linguistic abuse we call it a gaussian ap¬ 
proximation | Roughly speaking, this means taking the minimum value for the t* technical 
parameter, i.e. r* = 0. This formally relies on the assumption that integrating the system 
on the vector space defined by the tangent space at the critical point actually takes into 
account the relevant configurations giving the most important contribution to the functional 
integral. This was actually holding in the case of 0- 


Does the gaussian approximation hold true also for the CRM model? Figure [5] reveals 
that there is actually a regime in which it can do pretty well. Not surprisingly, it is a regime 
in which results are not that far away from the phase quenched approximation; we know 
that this is a regime in which the two different eigenvalues of the hessian at the critical point 
are quite close to each other and the problem appears all in all quite symmetric and not that 
far away from the regime of ([7]). Note that the value of the (rescaled) mass at which the 
solution provided by the gaussian approximation departs from the correct one varies with N. 


Next step was to leave the gaussian approximation aside and try to implement the full 
Aurora algorithm. There are a couple of issues one should be aware of: we need a solid 
estimate for r*; also, within a time length of order r* we have to make sure we have under 
good numerical control both the SA and the SD. The latter is the critical one, for which we 
have already made clear that a BVP formulation is the choice to go for. Our implementation 
was along the same lines of the code available at [3T] . All in all, our experience with the 
complete Aurora algorithm for the CRM model was at first somehow inconclusive: a clear- 
cut indication of values of t* at the same time safe and manageable was missing. We will 
come back to this observation later, in the framework of the other numerical approach that 
we chose to implement. 


4.3 A different numerical approach 

We now want to take advantage of the parametrization 

T tJW (n, t) 

(a few of the formulae we will need to implement our strategy were clearly stated in [3], 
while the strategy itself we will see at work in the following was first described in [32]). 

13 One should not confuse this with the gaussian approximation described in 3.1. In that case the action 
is approximated with its leading (gaussian) term and all the thimble analysis is performed consequently. In 
this (algorithmic) gaussian approximation one pretends that the thimble manifold we are interested in and 
the thimble associated to the gaussian approximation of the action sit on top of each other also away from 
the critical point. 
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Figure 5: Exact (solid red line), phase quenched (dashed blue line) and gaussian approxi¬ 
mation results for the condensate, at fixed Nf — 2, jl — 2. 


Basically one describes a generic point by locating it on the SA curve it lies on. This means 
providing n (the direction one is taking while leaving the critical point) and the time t at 
which one reaches <f> while integrating the SA equations. The Erst goal is now to rewrite the 
contribution to the partition function which is attached to one thimble. In full detail this 
reads 


Z ^ = / dzi A ... A dz n e s 

Jjv 

n 


s, IL n llj 

T / TT dyl det(C) e~ s — e~ iSl W / TT dy\ e“ e~ s ’ 

1_ JVc ■ , _. _ _ JVc 


(25) 


charts C c % 


charts C c i 


In (25) we have taken into account that Sj is constant on J a . Moreover, there could be more 


than one relevant chart and on a given chart we have to take into account the residual phase. 
For the sake of simplicity in notation, we now take a few shortcuts. First of all, we discard 
the overall phase e ~ tSl ; it will be easy to account for it when we come back to the actual 
computation of an observable. We also discard the fine detail of more than one chart, since 
in practice this is not a issucj^] Finally, we leave the residual phase aside, having in mind 
that we can take it into account a posteriori by reweighting. We thus write a new quantity, 
which is the one we will further manipulate in order to single out the contributions from the 


14 


We will be concerned with single ascents, for which we will revert to a different, smooth parametrization. 
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single ascents. We define^ 

/ n 

11^ (26) 

i =1 

Roughly speaking, this is the quantity that can have a probabilistic interpretation. The key 
point is now to write an expression for 1 


1 


/ n n. n 

JJ d n k 5 (|n| 2 — l) / dt JJ S (y, 

k= 1 i=1 


Vi(n,t)) 


(27) 


where {?/* (h, f)} are the coordinates of the field as expressed in the local (orthonormal) basis 
{[/x ( t )} parallel-transported along the SA defined by n until time t. The solution for ( t ) 
is in terms of (the module of) a determinant 


or 


A* (t) = 


( a(|nf-l) ^(|nf-l) 

St Sni 

S(yi-yi(n,t)) S(y 1 -y 1 (h,t)) 


det 


5t 


Sni 


S(yn-yn(n,t)) S(y„-y n (n,t)) 


An (t) = 


det 


^(Kl 2 -i) ^ 

5n n 

S{yi-yi(n,t)) 


5n n 


S(yu-yn(n,t)) 


St 

Sn± 

Sn 

( 0 

2n i 

2 Tin ^ 

Syi (n,t) 

Syi(h,t) 

Syi (n,t) 

St 

Sni 

Sn n 

Sy„(n,t) 

Sy n (n,t) 

Sy n {n,t ) 

\ St 

Sni 

Sn n ) 


(28) 


The first column of this determinant can be easily related to the gradient of the action. 
It turns out that to compute the generic matrix element we need to do the following: 


We need to evolve not only the field, but the entire basis by integrating i§. 

We construct the 2 n x n matrix V whose columns are the (t)}. 

We construct the 2 n x n matrix u whose columns are the vectors {nh )} i=1 n which are 
obtained from the {V"W (^) | by means of Gram-Schmidt orthonormalization procedure. 


The relation V = uE holds, with 


Eij — 


y(j) . u d) j > i 
0 j <i 


The entries of the determinant we are looking for are now given by 


f = E UKn t E. 


ik 


Sm 


(29) 


15 For the sake of notational simplicity we also omit the explicit indication that the integration is on the 
thimble, as it is easy to recognize we are assuming. Notice that (26) is the generalization of © of section 
3.1. 
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Not surprisingly, there is a lot of information in the (tremendous amount of) computations 
we have just sketched. In particular, if we now introduce the n x 2 n complex space projector 


P 


P — (inxn 



(30) 


then the n x n complex matrix U = Pu is unitary; this is precisely the matrix of eq. ([5]), 
whose determinant is the residual phase e luJ . 


The details of the previous computation of A* (t) are given in appendix [B] We now 
proceed to make use of the expression for the identity encoded in (27). Inserting it in (26) 
we get 


/ IO 

i= 1 

P n p n 

= I JJ dyi e~ SR A h (t) J JJdn fc <5 


n\ 


2=1 

n 


/ n 

dt n ^ t- h - yi 

i =i 


10 n n f O 

JJ d n k 5 (|?7| 2 - l) / dt / JJ dy t S (y t - y t (n, t)) A h ( t ) e 

k= 1 ^ ^ 2—1 

n p 

dn fc S (|?z | 2 - l) I dt A h ( t) 


D -S R {h,t) 


which has a possible interpretation in terms of 

/ n 

]Jdn k 6 


n\ 


1) 3 


(>) 

n 


(31) 


k =1 


i.e. there is a contribution to the partition function for each SA path 

+oo 


zg 1 = / diA ft («)e- s “ , ’ u) 


(32) 


Note that the procedure naturally defines a probability, i.e. that for a point reached at time 
t on the SA defined by n 

A - (A 

P, (t) = - (33) 

One can also naturally define the cumulative distribution function (it is manifestly non¬ 
decreasing, positive definite and has the correct normalization) 


F n\t) = ^P ) I dt'A h (t')e- S ^ 


Z, 


(34) 


Since we can easily invert this function numerically, we have a tool to ideally sample con¬ 
figurations on a single SA. Namely, we extract a random number £ e [0,1] and then get the 
point on the SA (rather, the time at which the point is reached) by t — F^ 1 (£). Actually 
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this is not that useful. The fact that we ascend all the way along a given SA in order to 
compute Z^ suggests that it is rather convenient to compute the entire contribution which 
is attached to that given ascent. On the other hand, the relative weight of a given SA (within 
the complete partition function Z^) is given by Z^ / Z^ a \ 


We now want to take advantage of the parametrization $ e X o (n, t ) in the compu¬ 
tation of an observable. In the following, we will assume we are in a case in which only one 
single thimble is relevant. This is not the general case, but for what we want to obtain it is 
not a limitation. In the cases in which more than one thimble contribute, we can address the 
problem using the same strategy described in 3.1 in the context of the quartic toy model: it 


will be easy for the reader to generalize eq. (14) and the discussion following it. With this 
caveat in mind, we first of all write 


(O) = 


I' dzi A ... A dz n O e 

“ Ua 


-s 


ZW 

f n;u dy* ° eiu e ~ SR 

z>icu p~Sr 


0)c 


(e iuj ) c 


/nr=i%e 

-Sr yy e j iave till now simply generalized eq. (13). 


... e 


where (... ) CT = Z (ff) 1 ./' ECU d ?/ 

We can go further by making use of the new parametrization we introducecf^l 

f Vn f dt A fi(t) e~ SR ( n,t ^ e l0J ( n ’ t ) 0(n,t ) 
f T>h f dt A n(t) e~ SR ( n ' t '> 

f Vn Z^ {z^ 1 J dt A n(t) e- SR( - n,t ' 1 0(n, t) 


<o> = 


f Vn Z^ i^zS^ 1 / dt A ft(t) e~ SR ( n ' t ') 
JVnZ^ (e‘“0) 4 


jVn Zf{e ■“) 


(35) 


Eq. (35) is in a sense a new average. Namely, the different directions n can now be regarded 


as the new degrees of freedom of the overall integral, the quantities to be measured are the 
(Oe“)- and (e^)^ (he. partial averages attached to single SA), and the weights are given 
by the Z^\ It is rather obvious that 


• The basic building block are complete ascents. This is good, since we can have their 
computation under good numerical control. In other words, sampling on the thimble 
is not a problem: we stay on the thimble by definition. 

• The way to importance sampling now appears tricky. This is easy to understand, since 
picking up a contribution means picking up a n, whose weight Z^ is not known a 
priori , but only after the SA path associated to h has been obtained. 

• The crudest approach one can think of is of course a uniform sampling of the h-space: 
this is a static, crude Monte Carlo, which can easily become inefficient (in particular 
for large systems). 

16 We denote Vn = nl~i <5 (jn| 2 — 1^ . 
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In the following we will just be satisfied with the last approach of static, crude Monte Carlo: 
this will be enough to show that we can reproduce the correct results for the model at hand 
(even in regions where the sign problem is quite severe) and this holds true taking into 
account the contribution of the single critical point we found. From this very basic approach 
there will be something to be learnt also with respect to the Aurora algorithm (and on the 
computation of density of states as well). We will finally report on a few speculations on 
smarter algorithms which we are trying to devise. 


4.4 Results for the CRM model 

Figure [6] displays the results we obtained from simulations performed in the static Monte 
Carlo approach we have just discussed. All these results come from the contribution of one 
single thimble. As we have already pointed out, one original motivation of ours turned out 
to be not relevant: we were ready for looking for dominance of one thimble in some asymp¬ 
totic regime (in the thermodynamic limit) and, in the region of parameters we studied, we 
actually found no other thimble but the trivial one. Also, complexifying the theory in the 
parametrization that was shown to be problematic for complex Langevin [12] did not result 
in any problem. Results are shown for AI = 1, 2, 3,4 and fixed Nj — 2, jx — 2. 



Figure 6: Exact (solid red line), phase quenched (dashed blue line) and thimble simulations 
results for the condensate, at fixed Nf — 2, jl — 2. 


It is interesting to regard the parametrization we have employed from another point of 
view. In the upper panel of figure 0 we plot the quantity A^(f) e s ^ n > t )/as a function 
of Sr along a given ascent (remember that on each ascent one single value of Sr is only read 
once). This is the real weight of the functional integral for the configurations which lie on 
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that given ascent: it can be thought of as a different way of looking at the density of states 
(namely, this is the contribution attached to a given ascent). 

In the lower panel of hgure[7]we plot the same quantity as a function of the flow time. A first 




-5 -4 -3 - 2-1 0 1 

flow time t 


Figure 7: Real weight of the functional integral for the configurations which lie on the SA 
defined by a particular n as a function of Sr (upper panel) and t (lower panel) for N = 2, 
m — 7, Nf — 2, fi — 2. 


remark is due for the long initial flat region. Consider eq. (jT]). When we want to compute 
the contribution from a single ascent (that is from a single n) we need an initial condition, 
i.e. an initial value to at which the asymptotic regime holds. In principle, the more back in 
time we take to, the better initial condition we prepare. There are of course accuracy issues 
one has to live with. The flat initial region reflects the fact that we do our best to ensure we 
stay on the thimble. On the other side, one would like to know till what value of the flow 
time the asymptotic regime holds to a reasonable confidence. We can think of more than 
one indicator for the latter condition, e.g. the gaussian approximation of the action is very 
close to the actual value of S, or the factor A ^(t) is very close to its gaussian approximation 
(see appendix |B.1[ ). We mark with a (red) star a value of flow time which can be assumed 
to be the boundary of the region we have just described. Now, an efficient dynamic Monte 
Carlo is supposed to sample configurations in regions where the weight is concentrated. In 
this sense we can say that the distance (in flow time) from the region around the maximum 
of A n(t) e ~ SR ^ n,t ^/and the flow time marked with the star is a reasonable indicator of 
the t* parameter of the Aurora algorithm. 


A comment is due concerning the residual phase: we encountered no problem in taking 
it into account by reweighting. As it was expected, it is a smooth function on the ascents, 
so that (Oe*")- and {e luJ ) fl can be safely computed and wild cancellations are never there 
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at any stage of our computations. The fact that the residual phase is smooth does not of 
course mean it has no net effect, as can be seen in figure [8j In the upper panel we show the 
effect of neglecting the contribution of the denominator of eq. (35): this amounts to com¬ 
puting Z^~ l f Vh Z^ (e lU} Q) f Pl In the lower panel we show yet another phase quenched 
computation, namely what we could term a residual phase quenched result. In this case 
we simply omit the contribution of the residual phase and compute Z^~ x J Vh Z^ (0) h . 
Both plots show that reweighting for the residual phase is essential to get the correct results; 
this happens to be the case in particular for low dimensions. 




Figure 8: The effect of not accounting for the residual phase. In the upper panel 
its contribution is not accounted in the denominator of eq. (35), i.e. we compute 
Z^~ l f Vh Z^ (e lu O)-. In the lower panel we omit the residual phase completely and 
compute Z (a ^~ l f Vh Z^ (0) h . 


17 Note that in this case we take the residual phase into account in the numerator. 
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We admittedly made use of the crudest possible application of the parametrization con¬ 
tained in eq. (J35|, i.e. a static, crude Monte Carlo sampling of the integral. It is static, i.e. it 
is not based on a stochastic process. It is crude, because it does not implement importance 
sampling^) In general, importance sampling computes an integral / = f dxf(x) as an av¬ 
erage I = J dx p(x) and all the point is being able to extract configurations distributed 
according to p(x). The natural importance sampling for (35) would try to sample the space 
of ascents (i.e. the different ft) according to the weights and we have already made 

the point that this is tricky. Crude Monte Carlo simply extracts n with constant probability. 
As it is well known, the more non-trivial the profile of Z^/Z^ is, the more inefficient one 
expects the crude Monte Carlo to be. This can be clearly seen at low masses, where the 
problem is less symmetric with respect to different choices of the ft: we have already made 
this point while commenting on the failure of the gaussian approximation for the Aurora 
algorithm. The interested reader is referred to appendix B.l to get more insight on this in 
the case of the gaussian action. Here we want to stress that the problem with crude Monte 
Carlo does not necessarily relate to the dimension N (of the matrices): going low enough 
in mass at any fixed N can already result in a difficult computatiorp^j As the mass gets 
lower, one indeed clearly sees larger error bars in figure |6j Needless to say, this is a region 
in which we had to collect many ascents: hundreds of thousands, actually, by definition all 
statistically independent. This is a huge numerical effort. On the other side, even accepting 
that at some point our crude Monte Carlo has to give up in front of a non-trivial profile of the 
weights Z^ / Z^ a \ it is reassuring to see that we nevertheless solved a non-trivial problem: 
figure [6] clearly shows that we were able to solve the problem also in regions in which the 
sign problem shows up as quite severe. 


All in all, we saw that implementing importance sampling is hard, since the relative 
weights Z^/Z^ are only known after the complete SA associated to n has been computed. 
This could sound a very pessimistic conclusion. On the o ther hand, the gaussian approxi¬ 
mation of the Z^ can be easily computed (see appendix B.l), which suggests the idea of 


making use of them to formulate proposals for the n. This is something we are currently 
investigating. 


5 Conclusions and prospects 


We discussed the solution of a simple toy model via thimble regularization. Quite interest¬ 
ingly, this model, which dates back to some thirty years ago and was proposed as a sort 
of benchmark for complex Langevin, was still missing a full solution in the context of the 
latter. In thimble regularization the solution is clear and can be implemented numerically 
by a number of simulation algorithms. 

18 We recall that the Monte Carlo methods we are mostly familiar with are dynamic Monte Carlo in which 
importance sampling is obtained via convergence to the equilibrium distribution of a stochastic process, 
which in virtually all the cases is a Markov chain. 

19 Of course going to higher values of N would result in extra computational effort, but what is really 
crucial is to see at each value of N what is the threshold in mass below which we have a too much non-trivial 
profile of the Z^ / Z^ a \ 
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We then investigated the Chiral Random Matrix model and showed that thimble regu¬ 
larization can successfully deal with the sign problem that the system displays. In the 
region of parameters we studied, a single thimble accounts for the results. We made use 
of a parametrization in terms of contributions attached to SA (which are the basic build¬ 
ing blocks to define the thimble). This was done by crude Monte Carlo, leaving open the 
problem of devising a smarter algorithm (importance sampling) to take advantage of the 
parametrization we made use of. This is the subject we are currently investigating in view 
of other applications. 
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A The Hessian for the CRM model 


We want to compute the hessian at the critical point a t j = b 7 j = a 7 j = f3ij = 0. We need the 
following second derivatives (the fields are complexified by setting e.g. a 7 j = a R + i oh etc.) 
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where use of the Cauchy-Riemann equations has been made (the other derivatives are 
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trivially related to these by Schwarz theorem, e.g. „ R £ R 

cfa mn^ a ij 

are given by 
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The hessian for N = 1 is (with the conventional choice of ordering: a R , b R , a R , f3 R , a 7 , b 1 , a 7 , (3 7 ) 
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As the second derivatives are manifestly diagonal with respect to the indices i,j,m,n, 
the hessian for a generic N is block-diagonal 


N 

tf(AO = 0tf(l) 


The model thus features a huge degeneracy of eigenvalues, as the spectrum of N = 1 is 
repeated N times. The spectrum for N = 1 features 4 positive eigenvalues and 4 eigenvalues 
opposite in sign (as expected by holomorphicity). We are interested in the positive part of 
the spectrum. An explicit computation shows that the distinct positive eigenvalues of 
are actually 2 and they are 


A± = 


2 m 2 


2Nf cosh (2p) ± y/2 \J 8m 4 — 8Nfm 2 + Nj + Nj cosh (4p) 


We note in passing that (here and in many other places) we could have written formulas 
in the complex notation of Takagi factorization theorem (see H), which we decided not to 
employ here and in all the paper to stick to a completely real notation. 


B Computing A^(t) 


We want to compute A a(t), which is defined in (27). The main point is that we are ascending 
along a given flow, and while doing that we are transporting along the flow also the basis 
vectors, i.e. we are integrating ([6]) as well. Near the critical point the ([7]) hold, in which the 
parametrization $ G J a (n,t) is manifest. Given a reference point t 0 1, the (JT|) can 














be regarded as initial conditions for the flow associated to n. Near a generic point, under 
infinitesimal variations of t and h, the variation of the point 5$ is given by 

n 

5$ = 5$ (n,t) = X^ W (t) 5c® 

1=1 

This is so because £$ is itself a vector belonging to the tangent space T§,J a . The (constant) 
coefficients 5c® can be worked out from the asymptotic form of <f> (2) near the critical point 
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from which it follows 


5c® = 5n.i + XiTiiSt 

Being a vector of we can write it as a decomposition on the (orthonormal) w-basis 
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and from this we have 
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Let us now consider the terms ^ appearing in (2), where * is either 2 or rij. For these 
we have 
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Now we make use of the explicit form of 5c^ k \ which gives = Aand = 5, 
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Jkj j 


from which one can easily derive the (29). 
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B.l A n(t) in the gaussian approximation 

We can compute A* ( t ) for the gaussian case (purely quadratic action), where the asymptotic 
form for the SA and parallel-transport equations is correct arbitrarily far away from the 
critical point. In that case the entries of A* (f) are (E = t nxn ) 


and the determinant is 
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For the gaussian action Sn(n, t) = A/? (by) + | Y^l=i ^k'nle 2Xkt , so that collecting every¬ 
thing we can write an expression for the weight itself 

From this expression it is easy to understand that the more the eigenvalues differ from each 
other (which in our case happens for low values of the mass parameter), the more various 
Zff* can differ. 


dt e 
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B.2 A useful consistency relation 

Since in the end we are performing quite a lot of computations, it is useful to have a consis¬ 
tency relation to be checked while ascending along the flow. The gradient of the action is 
yet another vector belonging to the tangent space T<^J a . Let us write the decomposition 


y(i) (*) 9 (i) 

i =1 


the coefficients (f l) can be found with the aid of the asymptotic form of the action 
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where we have used the symmetry of the hessian and the fact that Hv® = A iV^\ We 
have found that g w = so while integrating the flow equations, we can keep checked the 

norm 


i =1 


and make sure that it is small with respect to the size of the system. 
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